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The Hartree-Fock-Bogoliubov equation for the ground states of even-even atomic nuclei is solved using the 
canonical representation in the coordinate space for zero range interactions like the Skyrme force. The gradient 
method is improved for faster convergence to the solutions under constraint of orthogonality between canonical 
orbitals. Necessity of the cut-off of the pairing interaction is shown even when the number of the canonical orbitals 
are restricted. A repulsive dependence of the interaction on the pairing density is introduced as an implementation 
of the cut-off which leaves the HFB super matrix state-independent. 



Introduction 

In neutron-rich nuclei, the pairing correlation significantly 
involves the continuum single-particle states. This makes 
the HF+BCS approximation inadequate due to the un- 
localization of the neutron density distribution and de- 
mands one to solve the Hartree-Fock-Bogoliubov (HFB) 
equation without approximations. The solution, in the 



coordinate space was first formulated in Ref. u using 
the quasi-particle states and performed for spherically 
symmetric states. However, its application to deformed 
states is difficult because there are quite a large number 
of quasiparticle states even for a moderate size of the 
normalization box (i.e., the cavity to confine the nucle- 
ons to discretize the single-particle states). Every HFB 
solution has an equivalent expression of BCS variational 
form. The single-particle states to construct the BCS 
type wavefunctions are called the HFB canonical basis 
or sometimes the natural orbitals. This expression was 
used_to solve the HFB for spherical states originally in 



Ref.B 



Although spherical solutions can be obtained easily with 
present computers (for zero-range forces), deformed saht 
tions are still difficult to obtain. The two-basis methodlmcP 
is the only one implemented so far for neutron-rich nu- 
clei. Some recent developments like Reflf are also in 
progress. 

We have applied the canonical-basis method to deformed 
states in a three dimensional cubic mesh-representation 
with density dependent delta interactions.^ It has turned 
out to be a very efficient alternative approach to obtain 
the solutions. The origin of its effectiveness is that the 
number of necessary single-particle basis states to de- 
scribe the ground state of a nucleus is proportional to 
the volume of the nucleus in the canonical-basis method 
while it commensurates with the volume of the normal- 
ization box in conventional methods. The difference of 
the number of the basis states amounts to a factor of 10 1 
- 10 3 . 



In this paper we discuss the canonical-basis formulation 
of the HFB, the method to obtain the canonical-basis so- 
lutions, faster gradient-method paths than a naive imag- 
inary-time evolution, the necessity of the cut-off of the 
pairing interaction, and an implementation of the cut- 
off in terms of an interaction dependent on the pairing 
density. 



HFB in the canonical representation 

To begin with, let us formulate the HF and the HFB in 
the coordinate-space representation in order to elucidate 
a difficulty of the HFB and to suggest its possible solution 
in terms of the canonical-basis representation. For the 
sake of simplicity, we consider only one kind of nucleons 
and designate the number of nucleons by N in Eqs. (0)- 
(0), which are in this section. The z-component of the 
spin of a nucleon is represented by s (= ±i). 

In the HF, one should minimize (9\H for single Slater- 
determinant states, 



N 



m = n a !i°>' 

at = ^2 J drij]i{r,s) (J{r,s), 



(1) 
(2) 



by varying {V't}t=i,"-,JV under orthonormality conditions 
{ipi\ipj) = Sij. The operator aj creates a nucleon with a 
wavefunction ipi(r, s). The distribution function of the 
density of the nucleons is related to the wavefunctions as 



N 



p(r)=^(*|a(f, S )at(f, S )|*) = ^^|^(f )S )| 2 .(3) 
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In the HFB, the solution takes the following form, 



I*) 



n>io>, 

^2 / dr{<f>i(r,s) a(r,s) 

s J 

+(fi(r,s) a\f,s)} , 



(4) 



(5) 



where bi is the annihilation operator of a negative-energy 
Bogoliubov quasi-particle with amplitudes tfii (r, s) for pres- 
ence and (pi(r,s) for absence. / is the number of the 
basis states of the employed representation. For a threa- 
dimensional Cartesian mesh (3D-mesh) representation^ 
it is the number of the mesh points (times four when 
spin-orbit interactions are included) and typically 10 4 - 
10 5 . One should vary {ipi, ^i}i=i,...,i under orthonormal- 
ity conditions 



s J 

% (i<*<i<i), (6) 



and a constraint on the expectation value of the number 
of nuclcons, 



(*\N\V) = [ dr |^(r,s)| 2 =iV, 

1=1 s ^ 

iV = y dra(r, s)^ (r, s). 



(7) 
(8) 



where a] and a| create a nucleon with wavefunction tpi (r, s) 
and i/^r-s), respectively, which ape called as the canon- 
ical basistP or the natural orbitalstJ of the HFB vacuum 
\ty). One must use K = i/ for the exact equivalence 
between Eqs. (|4|) and ( |To| ) in general cases. When \^>) is 
a time-reversal invariant state, which we assume in this 
paper, ipi and ip\ are the time-reversal state of each other. 
In this case, only one wavefunction of each time reversal 
pair should be counted as independent variables of the 
variational procedure. 

To obtain solutions in the canonical-basis framework, one 
should vary {ipi, Uj, «i}j=i,...,.K' under three kinds of con- 
straints, i.e, the orthonormality conditions, 



s J 

= 5 l3 (l<i<j<K), 
fixed expectation value of the number of nucleons, 

K 

i=i 



(12) 



(13) 



and the normalization of the u-v factors + vf = 1. 



The nucleon density is expressed as 



K 



(14) 



The nucleon density for state (^) is expressed as 
i 

i—l s 

The essential difference between the HF and the HFB 
is that one has to consider only N ~ 10 2 wavefunctions 
in the former while one has to treat explicitly as many 
single-particle wavefunctions as the number of the basis 
in the latter. 

Owing to the Bloch-Messiah theorem,^ the state (||) can 
be expressed (for the ground states of even-even nuclei) 
as, 

K 

i*) = n( u * +w *°N)i°>' ( io ) 

i=l 

at = ^ dr'4>i(r,s)rf(r,s), (11) 



Rcinhard et al. regarded that the advantage of the repre- 
sentation ([h]) over (|J) is that one has to consider only a 
single set of wavefunctions {ipi}i=i,.-.,i unlike a double set 
{<j)i, fi}i=±i,---.±i/2" > However, we expect much greater 
benefit from the canonical-basis representation. Namely, 
i may be truncated as i < K — 0{N) «C \l to a very 
good approximation. It is because ipi appearing on the 
right-hand side of Eq. ([l4]) must be a localized function 
as p{r) on the left-hand side, while the orthogonality ( |l2] ) 
does not allow many low-energy wavefunctions to exist in 
the vicinity of the nucleus. For 3D-mesh representations, 
/ is proportional to the volume of the cavity while K is 
proportional to the volume of the nucleus. The latter is 
lf^-lO 3 times as small as the former. 

Incidentally, the situation is quite different in the quasi- 
particle formalism. On the one hand, the localization of 
the density demands only the localization of ipi through 
Eq. @ while does not have to be localized in general. 
On the other hand, the orthogonality condition ^ in- 
volves both (fi and (pi. This discrepancy enables many 
quasiparticle states having similar tpi to be orthogonal to 
each other by differing their tpi. 
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Mean fields for zero-range interactions 

Let us present the effective Hamiltonian employed in this 
paper. We adopt a density-dependent zero-range interac- 
tion. Zero-rangeness makes the mean-field potentials lo- 
cal, which is an essential advantage for coordinate-space 
solutions. On the other hand, the omission of momen- 
tum dependences are merely for the sake of simplicity 
and there will not be essential differences i» the formula- 
tion if we use the full-form Skyrme force.EP Our force is 
expressed using the parameterization of the Skyrme force 
as, 



v(ri, si; fa, s 2 ) 



(15) 



where fj and Si are the position vector and the spin state 
of the two interacting nucleons, i = 1,2. Dependence 
on the isospin state is redundant because of the zero- 
rangeness. We adopt t = -983.4 MeV fm 3 , t 3 = 13106 
MeV fm 3+3ct , and a = 0.98 when the -force is used to con- 
struct the mean- field (HF) potential.^ We use different 
strengths to make the pairing potential. Wei-express the 
pairing force in the parameterization of Ref.tiP as 



1-P< 



l-^Ww 2 ). (16) 

Pc J 



P a is the exchange operator of spin variables. ^(1 — P a ) 
is a projector into spin-singlet states so that the interac- 
tion acts only between like nucleons. We use p c — 0-32 
fm~ 3 (to roughly vanish the volume-changing effecttiP) 



and v p = -440 MeV fm 3 . 



When one considers both protons and neutrons, the state 
of the nucleus is expressed as a product of two BCS type 
states (|l^) for the protons and the neutrons: 



spin, the wavefunctions ipi(r, s) can be factorized into a 
product of a spin wavefunction and a real function of the 
position, which we write ipi(r) in the following. It holds 



With the interactions ([l5]) and (|l6|), the total energy for 
state dl7]) is written as, 



(18) 



E = = J H{r)dr, 



(19) 



where m is the average of the proton and the neutron 
masses divided by 1 — 1/A for the correction of the center- 
of-mass motion. H(r) is called the Hamiltonian density 
while function of position r in the right-hand side are 



K 



r(r) = g^^v 2 \Vipi(r)\ 2 : kinetic energy density, 
»=i 

K 

P(r) = g^2v 2 \^(r)\ 2 : density, 
z=l 
K 

P~(r) = 9^2,UiVi\ipi(r)\ 2 : pairing density, 



where 5 = 4, which is a factor to account for the situation 
that a wavefunction ^ takes care of four nucleons for the 
spin and isospin degeneracy. 

The mean-field potential V and the pairing potential V 
are defined as 



V = = -hp + —rrrhp - -r^-p, 
op 4 16 8p c 



(20) 
(21) 



i*>=nn 



Ui + v t at q aj, q 



|0>, 



(17) 



q=p i=l 



where q distinguishes between protons (p) and neutrons 
(n). a\ creates a proton having a wavefunction ipi !P (r, s) 

while a\ creates a neutron with a wavefunction ipi <n (f, s). 
The product form is due to the pairing interaction ( fig ) 
acting only between like nucleons. 

For the sake of simplicity, we treat N=Z nuclei with- 
out Coulomb interaction in this paper. In this case, the 
wavefunctions are the same between protons and neu- 
trons, i.e., ipi, p (f,a) = ip itU (r,s), ip^ p (r,s) = ipi >n (r,a). 
Moreover, because the potentials are independent of the 



in which V, V , TL, p, and p are local functions of r while 
to, t 3 , a, v p , and p c are constants. 



The mean-field and the pairing Hamiltonians are 



h 
h 



^V 2 
2m 



V. 



V. 



(22) 
(23) 



The quasiparticle states are the eigenvectors of the so- 
called HFB super matrix composed of h and h: 



(24) 
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This is just an eigenvalue problem of a hermitian (be- 
cause of the time-reversal symmetry) matrix. The canon- 
ical orbitals are also determined by h and h but in a more 
complex way as described in the next section. 



Gradient method for canonical-basis HFB 

In this section we describe a procedure to obtain the 
canonical-basis solution of the HFB equation directory, 
not by way of quasi-particle states. Instead of mini- 
mizing E = {^\H\^} with |*) given by Eq. (0) under 
constraints of Eqs. ( p"2| ) and ([l3|), one may introduce a 
Routhian R, 

K 

R = E-e F -gY,v} 

i=i 

K K 
<=1 2=1 

and minimize it without constraints, ep is probably the 
most familiar Lagrange multiplier, whose physical mean- 
ing is the Fermi level. In the definition K 2 Lagrange 
multipliers \j obeying hermiticity, 

Kj = A*„ (26) 



following equation: 

1 SR „ , , ^ , , 

% 3=1 
K K „ 

- EE^{^w-w = °> ( 29 ) 

J= i k=i Vi 

where 

Hi = vfh + UiVih. (30) 

One can regard Hi as a state-dependent single-particle 
Hamiltonian. This dependence on states makes the or- 
thogonalization conditions essential to the method. For 
HF, the orthogonalization conditions are easily fulfilled 
because tpi are eigenstates of the same hermite opera- 
tor h and thus are orthogonal between themselves at the 
solution : {ipj\ipi) ■ (e,- — e,) = 0. The orthogonaliza- 
tion procedure is needed only because states satisfying 
the orthogonality are unstable for decaying into Pauli- 
forbidden configurations. On the other hand, for the 
canonical-basis HFB method, the orthogonalization is es- 
sential because the single-particle Hamiltonians Hi differs 
from state to state. Therefore, the determination of the 
explicit functional form of Xij is the most important part 
of the method. Reinhard et al. have proposed 



are introduced instead of \K{K + 1) independent mul- 
tipliers. This hermitization of A is adopted in order to 
make R real so that two conditions, SR/Sipi = and 
SR/Sip* = 0, become equivalent and thus one has to con- 
sider only one of them. Note that 5{j is subtracted from 
(ijji\ipj}, in contrast to Ref. EP, which is in order to treat 
Xij not as constants like eF but as functionals of the wave- 
functions. 



Ay = {Hi+H^i). (31) 

Let us reason on which grounds the above definition can 
be deduced. Understanding these grounds is indispens- 
able in order to modify the definition later for faster con- 
vergences. From the requirement that Eq. (|29| ) must hold 
at the solution (where (ipi\ipj) = one can deduce, 



Using notations, 



Xij = (ipj\Hi\ipi) at the solution. 



(32) 



(V>iW^>, A i = -(ip i \h\ip i ), 



(27) 



the stationary conditions of R result in two kinds of equa- 
tions. One is dR/dvi = 0, which concerns the occupation 
amplitudes Vi and is fulfilled by 



2 11 £i ~ £F 

vj = — ± — 

2 2^ 



(28) 



The minimum of R for the variations of Ui and Vi is ob- 
tained (when Aj > 0) by taking the minus sign for the 
double sign in the right-hand side of E q. (|28|) and choos- 
ing the same sign as Vi for Ui = ±yl — vf. The sta- 
tionary condition for a wavefunction ipi(f,s) gives the 



Eqs. (^Tj) and (32) are equivalent at the solution because 
X^ is defined to be hermite by Eq. (26). Since this her- 



miticity must hold at any points to ensure the equality 
between the number of constraints and the number of 
independent multipliers, one should not adopt Eq. ( |3^ ) 
but Eq. d|). 

One can utilize the gradient method to obtain the HFB 
solutions in the canonical-basis formalism. The most 
naive implementation of the gradient method agrees with 
the imaginary-time evolution method in its first order ap- 
proximation of the size of the imaginary-time step At: 



ipi -»• A 



1 a SR 
-At 

9 Sip* 



(33) 
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We have developed a 3D-mesh canonical-basis HFB pro- 
gram from scratch according to the above formulation. 
We take an example of our calculations using the pro- 
gram for the ground state of 40 Ca. The wavefunctions 
are expressed with 39 x 39 x 39 mesh points with mesh 
spacing of a=0.8 fm. We employed the 17-point finite- 
difference approximation to the Laplacian. Note that 
the requirement of precision is higher for HFB than for 
HF because one has to treat larger momentum compo- 
nents than the Fermi momentum in HFB. The vanishing 
boundary conditions are imposed on the boundary (the 
Oth and the 40th mesh points) and the wavefunctions are 
anti-symmetrically reflected in the boundary to apply the 
finite-difference formula. We considered K = 20 canoni- 
cal basis, which can contain Kg = 80 (=2 x A, A = 40) 
nucleons. 



For the imaginary time step size At, it must hof 



2 h 

^ ~ : J- max — 3 



dEi 



We took Ar = 1/T n 



2m \a 



(34) 



In Fig. H, the error of the second equality of Eqs. (|2^) ne- 
glecting the error of orthogonality, i.e., maxi = i i ...,x|HiV ; i — 
Yl - Xijipj\ is plotted with a solid line versus the number 
of evolution steps. The corresponding quantity for HF, 
max i=lj ... jj 4/4|/i'(/>i — {ipi\h\ipi)ipi\ , is also plotted with 
a dash line. The figure demonstrates that one can in- 
deed obtain HFB solutions with the natural-orbital HFB 
method in the 3D-mesh representation. We obtained sim- 
ilar convergence curves for the error of the orthogonality 
and for the inconsistency between the potential and the 
densities. 

The speed of the convergence is, however, about ten times 
as slow as the HF case. This is the subject of the next 
section. 



Figure [I] 



Incidentally, the period of 25 steps may be too frequent 
because the effect seems to saturate at periods around 
100. We adopt the period of 25 in most calculations, 
however, because the increase in the computation time is 
only a several percent of the total time with this period. 



Acceleration of the gradient method 

We show the origin of the slow convergence and present 
a solution of the difficulty in this section. 

Stccpcst-dcsccnt paths, which the gradient method draws, 
depend on the choice of the independent variables. For 
example, Eq. ( |33| ) is obtained when one uses (Re ipi, Im ipi) 
as independent variables to define the gradient vector and 
then express it in coordinates (ipi,ip*). If one uses scale- 

— 1/2 

transformed wavefunctions Xi = a i w%i where on is a 
scaling factor, a gradient step becomes, 



1 A SR 

9 Sxt 



which is equivalent to 



1 a 

Ipi -> ipi a i^ T TT7- 

9 oijj* 



(35) 



(36) 



The change from Eq. ( p3|) to Eq. (36) is equivalent to 
multiplying a,; to the single-particle Hamiltonian TCi in 
Eqs. A 



When one parameterizes the scaling factor as on = v\ 
the modified single-particle Hamiltonian becomes 



-2i/ 



,.2-2i/ 



l-2i/„ 



h + v , 

vfh + ViUih (v = 0), 

Vih + Uih (v = i), 

h + %h (y = \). 



(37) 



Fig. 1. Convergence to the HF and HFB solutions for 
40 Ca. 

We adopted an additional procedure which is not indis- 
pensable to obtain the solutions but effective to make 
the convergence more robust and somewhat quicker: Af- 
ter every 25th gradient-method steps, we orthogonalize 
{ipi} with the Gram-Schmidt algorithm in the ascending 
order of e,, defined in Eqs. (p7|), and then diagonalize the 
super matrix of the HFB Hamiltonian (^J) by expand- 
ing the quasi-particle wavefunctions (<fri,<Pi) (1 < i < K 
and their negative-energy partners) in a 2_ftT-dimensional 
basis {ipi} © {ipi}, and finally transform the resulting 
quasi-particle wavefunctions to canonical orbitals and oc- 
cupation amplitudes {ipi,Vi} to renew them. 



When v = (i.e., 014 — 1), to which the imaginary-time 
evolution ( |33|) corresponds, the single-particle Hamilto- 
nian aiHi {=Hi) can be very small for canonical orbitals 
whose €i is much higher than the Fermi level (i.e., ej — 
ep 3> A^. This smallness makes the changes of such 
orbitals very slow. On the other hand, for v = 1, the po- 
tential can be very deep for such high-lying orbitals due 
to the factor m / Vi in front of ft,. In this case, however, the 
gradient step may be numerically dangerous. We usually 
use ^ < v < 1, which provides a fast and numerically 
stable method of solution. 

When one introduces the acceleration factors (i.e., on > 
1), the multipliers Ay should be modified from Eq. ( |3l| ) 
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by the following reason. In the computation of a gradient 
vector given by the first of Eqs. (^9|), the last term takes 
much more computing time than the first two terms due 
to 8Xjk/6ip*. One can forget the last term if the orthog- 
onality relations ( |l2| ) are fulfilled along the path of the 
steepest descent. Let's suppose that the relations are sat- 
isfied before a gradient-method step is taken and require 
that they are conserved to the first order in At after the 
step, i.e., 



(VHVj) = Sij + O ((At) 2 ) if (frWj) = Sa, 



(38) 



of deformation 7. The last two quantities are determined 
from the mass quadrupole moments 2$ 

The dot curves were obtained without accelerations, i.e., 
with cii = 1 or v — in Eq. (37), while the solid ones were 
obtained with the acceleration method with cc,- = 1/ v- L 
or v = |. One can see that the convergences of these 
quantities become by far faster by using the acceleration 
method. This result demonstrates that canonical-basis 
HFB can be solved without very heavy numerical com- 
putations. 



with 



Cut-off of the pairing interaction for canonical-basis HFB 



ip'i = ipi 



(39) 



Finally let us discuss on the cut-off schemes of the par- 
ing interaction in relation to the canonical-basis HFB 
method. 



Substituting Eq. (|39j) into Eq. 
miticity (Eq) result in 



and requiring the her- 



1 



a,: + a 



■(tpj\ (ctiHi + OLjHj) \ipi). 



(40) 



Delta function forces without cut-off leads ta a divergence 
of the strength of the pairing correlation.E 3 } In order to 
circumvent the divergence, in the conventional method 
of solution, one usually takes only quasiparticles whose 
excitation energy e qp is lower than some cut-.off energy 
parameter e cu t to construct the ground state.El Namely 
Eq. (0) is modified to 



This form of Ay fulfills the requirement that it should 
agree with the expression (|3^) at the solution as the naive 
form of Eq. ([5l]) does. Two forms differ, however, before 
reaching the solution if one assumes an 7^ aj in general. 
Therefore, in order to conserve the orthogonality to van- 
ish the last term in Eq. (^9|) one must not use Eq. ( pl| ) 
but Eq. (|4(]) . We have indeed suffered from large errors of 
orthogonality by using the naive form ( |3l| ) . On the other 
hand, by using the correct form (|40|) , we have observed 
not only that the error does not grow during the evolu- 
tion but also that the error decreases without performing 
explicit orthogonalization procedures periodically during 
the evolution. This decrease should originate in the sec- 
ond order terms in At in Eq. (^8|), whose effects we did 
not consider. 



Figure 



Fig. 2. Convergence to the HFB solution for 3 S. 

We compare the results of calculations between v = 
and v = i in Fig. g. The wavefunctions are expressed 
with 19 x 19 x 19 mesh points with mesh spacing of a=1.0 
fm. We considered K = 16 canonical basis, which can 
contain Kg — 64 (=2 x A) nucleons. The figure shows the 
convergence history to a HFB solution. Four quantities 
are plotted as functions of the number of gradient steps. 
They are, from the top to the bottom, the total energy E, 
the pairing gap A (averaged with weight UiVi), the size 
of quadrupole deformation (3, and the size of triaxiality 



,qp\ 



6i|0) 



(41) 



where 9 is the step function. 

In the canonical-basis method, the restriction on the num- 
ber of canonical orbitals may prevent the divergence with- 
out introducing explicitly a cut-off energy. We have ex- 
amined this idea by performing numerical calculations 
for various situations. Then, we noticed that observables 
sometimes jumps suddenly in the course of long-time evo- 
lution. 



Figure | 



Fig. 3. Sudden changes in three quantities during the 
course of a gradient-method evolution. 

An example is shown in Fig. |^. The calculation was done 
for a nucleus 32 S on a 19x19x19 cubic mesh with a mesh 
spacing a=l fm. We considered K = 20 canonical or- 
bitals, to which we gave harmonic oscillator wavefunc- 
tions at the beginning. An 11-point formula was em- 
ployed for the Laplacian. The acceleration parameter in 
Eq. (37) is taken as v = 0.7. At first we suspected that 
these jumps were due to the acceleration method. How- 
ever, with smaller u, we still observed jumps; only they 
come later. 
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We investigated the origin of the jumps and found that 
each sudden change was due to a shrinkage of a high- 
lying (i.e., having large ej) canonical orbital to a mesh 
point. This shrinkage can decrease the total energy of 
the nucleus for the following reason: The contribution 
of a canonical orbital to the kinetic energy density is 
proportional to its BCS occupation probability vf , while 
pairing density is proportional to UjUj. Because Vi <C 1 
and m = 1 for high- lying orbitals, it holds that mvi ^> vf. 
Therefore the increase in the kinetic energy due to the 
shrinkage is easily compensated by the gain of the pairing 
correlation energy at the shrunken point. 

The observed jumps indicate that, in most cases (or maybe 
all the cases), there are no potential barriers between 
such physically meaningless solutions and the physically 
reasonable one. 

Incidentally, we confirmed that the lack of potential bar- 
riers was not due to the discrete approximation of the 
Laplacian in the kinetic energy term: The approxima- 
tion was based on the Lagrange polynomial interpolation, 
which tend to underestimate the expectation value for 
high momentum components. One might suspect that 
more accurate treatments of the kinetic energy could 
restore a barrier between the physical and unphysical 
solutions. However, we observed similar jumps by us- 
ing the Fourier transformation with periodic boundary 
conditions JiiP which gives the exact result up to n/a. 

We have decided that it is necessary to introduce cut-off 
for the reliability of the method. 

As the cut-off scheme, we first employed cut-off factors 
which are dependent on orbitals. n In the BCS approxi- 
mation, a smooth cut-off methods' is often utilized, in 
which the interaction is modified as 



with 



fi = CX P 



k? 



1.2fm, 



(45) 



We made the dependence to be on the kinetic energy 
(oc kf), not on the mean-field energy as in Eq. ([43|), 
to avoid a highly complicated expression of the gradient 
for the latter case. (In BCS, such complications are just 
neglected.) 

The result was successful to prevent the shrinkages to 
points. However, this cut-off scheme has an disadvantage 
that the HFB super matrix in Eq. (|4|) cannot be defined. 
It follows that we do not have well-defined quasiparti- 
cle states and cannot utilize them to express the HFB 
ground state. This drawback is rather serious because 
quasiparticles are useful to improve the precision of HFB 
solutions obtained by the canonical-basis method and to 
accelerate the convergence further, as described in the 
fourth section. 

As an alternative method, we have introduced a repulsive 
pairing-density dependence to the pairing force, Eq. (|l6|), 
in addition to the usual density dependence: 



v p (fi,si;f2,s 2 ) = v 




5 (fi - fa) • 



(46) 



A set of reasonable values of the parameters are v p 
-440 McV fm 3 , p c = 0.32 fm" 3 , and p c = 0.3 fm~ 3 . 



•^pair 



(\ / \ With forces ( plf ) and (46), the expectation value of the 

fi a\a\ I I fj aja , (42) energy for N — Z systems is expressed as a space integral 
i ) \ j I of a Hamiltonian density: 



fi = ffa), 



(43) 



where /(e) is a function of single-particle energy e. The 
function takes on ~ 1 well below a chosen cut-off energy 
and smoothly becomes zero above it. In analogy to the 
smooth cut-off method, we modified the pairing density 
as 



K 



K 



P = UiVi \^i\ 2 —> fi UiV t \ipi 



(44) 
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i=l 



+ 




(47) 



The mean-field potential V remains the same as Eq. (|2 
while the pairing potential V has an additional term 



V = 



dn 
"dp 7 



(48) 



With this new type of force, the shrinkage problem is 
completely removed. 



Figure [| 



plateau before the divergence sets in 



Fig. 4. Effect of changing the pairing-density cut-off pa- 
rameter p c (fm~ 3 ) on the pairing gap averaged with 
weight factor of UiVi for 32 S. 

In Fig. U we show the dependence of the pairing gap 
(averaged with weight UiVi) on the parameter p c , which 
controls the pairing density dependence. The values are 
taken after 5,000 gradient steps. The set up of the calcu- 
lations are the same as in Fig. [| except for p c . One can 
see that the pairing gap has reasonable values with p c < 1 
fm~ 3 . This is a good news because the pairing-density- 
dependent term, ~(p/p c ) 2 , can be small for the values 
of p which one finds in physical solutions, in contrast to 
the usual density-dependent term, —p/p ci which cancels 
roughly 50% of the density-independent term inside nu- 
cleus. This situation is illustrated in Fig. [|. The plotted 
quantity is the pairing Hamiltonian density 7i, which is 
the term in the second line of Eq. (47). The introduction 
of the new term demands only little change of the other 
parameters of the force if one adopts p c = 0.3 fm -3 . 
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Fig. 6. Evolution history of e», i.e., the expecta- 
tion values of the mean-field Hamiltonian for canon- 
ical orbitals. The pairing-density dependent term is 
switched on in intervals I and III (p c — 0.3 fm -3 ). It 
is suppressed to be very weak in interval II (p c = 3 
fm- 3 ). 



Conclusions 

We have developed a method to obtain canonical-basis 
HFB solutions in a coordinate-space three-dimensional 
(3D) Cartesian mesh representation. The features of our 
method are summarized as follows. 

i) It is not for spherical but for deformed nuclei and thus 
it can treat both deformation and continuum pairing si- 
multaneously. 

ii) It is not based on the oscillator-basis expansion but 
described in the 3D Cartesian mesh representation and 
thus can treat e.g. deformed halo-like orbitals. 



Fig. 5. Dependence of the pairing Hamiltonian density TL 
on the pairing density p for two values of a parameter 
p c =0.3 fm~ 3 (solid curve) and p c —oo (dot curve). The 
vertical scale is arbitrary. These curves applies when 

P < Pc- 

We show an example of the time evolution of in Fig. pi 
The set up of the calculation is the same as in Figs. [| 
and |] except that a 7-point approximation is used for the 
Laplacian to favor the emergence of unphysical solutions 
and make this calculation a very severe test of the cut- 
off scheme. The abscissa is the number of the gradient 
steps while the ordinate is the expectation value q of the 
mean-field Hamiltonian for all the K canonical orbitals. 
The pairing density dependence parameter p c is set at the 
standard value of 0.3 fm -3 before step 3,000 (interval I) 
and after step 6,000 (interval III). Between steps 3,000 
and 6,000 (interval II) p c is increased temporally to 3 
fm~ 3 , which results in a too weak dependence to prevent 
the divergence. One can see that an unphysical solution 
emerges without the pairing-density dependence term (in 
interval II) and that the unphysical solution is suppressed 
(in interval I) and is quickly restored to a physical one 
(in interval III) with the presence of the term. 

Before ending the section, let us mention that our next 
subject concerning the cut-off is whether repulsive mo- 
mentum dependences of the Skyrme force help to make 
the cut-off scheme more natural because it has been known 
that such momentum dependences make well-developed 



iii) There is a strong reason to believe that the necessary 
number of canonical orbitals is much smaller than the 
number of single-particle basis. 

iv) In order to perform variations under constraint of or- 
thogonality between the wavefunctions of the canonical 
orbitals, Lagrange multipliers werejjntroduced as func- 
tionals of the wavefunctions in Ref.B . We have clarified 
the necessary conditions for the form of these Lagrange 
multiplier functionals. We have modified the function- 
als appropriately so that they conserve the orthogonality 
during the course of the accelerated evolutions explained 
in the next item. 

v) We have found that the convergence to HFB solutions 
is very slow when one employs a naive gradient method. 
On the other hand, the convergence is quite rapid for 
Hartree-Fock solutions, which neglects the pairing corre- 
lation. We have investigated the origin of the slow con- 
vergence and found that the time scale of the gradient 
evolution is different from one canonical orbital to an- 
other depending on their BCS occupation amplitudes Wj. 
The difference ranges over many orders of magnitude. 
We have introduced an orbital-dependent acceleration 
method of the gradient evolutions and could overcome 
the difficulty of the slow convergence. 

vi) We have examined the effects of the cut-off of the 
pairing interaction. The 3D mesh mean-field methods 
have a practical use only for zero-range interactions like 



the Skyrme force presently. One needs a cut-off for zero- 
range pairing forces, without. which the pairing correla- 
tion energy diverges to — ootB We have found that zero- 
range forces need cut-off even when the number of canon- 
ical orbitals are finite. The divergence occurs through 
shrinking of high-energy canonical orbital(s) to (a) mesh 
point (s). 

vii) To suppress the divergence, we have introduced a 
pairing-density dependent interaction as a better choice 
than orbital-dependent cut-off factors. 

We believe that, by choosing a faster gradient path with 
the acceleration method and adopting the cut-off scheme 
in terms of the pairing-density dependence, the canonical- 
basis HFB method is now fully understood and has the 
potential to become the standard method to treat neutron- 
rich nuclei ixL.HFB. The details of this work will be pub- 
lished soon. I13 
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